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Abstract : 

The  passive  scalars  in  the  decaying  compressible  turbulence  with  the  initial  Reynolds  number 
(defined  by  Taylor  scale  and  RMS  velocity)  Re=72,  the  initial  turbulent  Mach  numbers  (defined  by 
RMS  velocity  and  mean  sound  speed)12  ’1  Mt=0.2-0.9,  and  the  Schmidt  numbers  of  passive  scalar 
Sc=2-10  are  numerically  simulated  by  using  a  7th  order  upwind  difference  scheme  and  8th  order  group 
velocity  control  scheme.  The  computed  results  are  validated  with  different  numerical  methods  and 

different  mesh  sizes.  The  Batchelor  scaling  with  k~l  range  is  found  in  scalar  spectra.  The  passive 
scalar  spectra  decay  faster  with  increasing  the  turbulent  Mach  number.  The  sheet-like  structures  of  the 
passive  scalar  gradient  are  also  founded.  The  extended  self-similarity  (ESS)  is  found  in  the  passive 
scalar  of  compressible  turbulence. 
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1.  Introduction 

Direct  Numerical  Simulation  (DNS)  becomes  an  important  tool  in  recent  research  of  turbulence. 
DNS  of  compressible  turbulence  is  more  difficult  than  that  of  the  incompressible  turbulence.  When  the 
turbulent  Mach  number  is  greater  than  0.3  the  shocklets  may  appear  in  the  compressible  turbulent  flow 
fields.  The  reason  and  mechanism  of  shocklets  existence  is  not  clear  yet.  The  turbulent  Mach  number 
in  DNS  can  not  be  very  high  with  the  present  existing  numerical  methods  and  computer  resource.  For 
the  problem  of  compressible  isotropic  turbulence  with  the  initial  Reynolds  number  Re=  72  (based  on 
Taylor  scale  and  RMS  velocity),  as  the  initial  turbulent  Mach  number  is  less  than  0.5,  the  tenth  order 
accurate  symmetrical  compact  difference  method  in  Ref  [2]  can  work,  but  as  the  turbulent  Mach 
number  is  greater  than  0.5,  this  method  is  not  valid. 

The  same  problems  but  with  much  higher  turbulent  Mach  numbers  Mt=0.7  and  0.95  are  solved  by 
authors  of  this  paper  with  high-order  upwind  difference  scheme  and  GVC8  scheme1’  41. 

There  are  many  problems  related  with  the  scalar  flux  in  turbulence,  such  as  the  pollutant  density 
in  the  air,  chemical  or  biological  species  concentration  and  salinity  in  the  ocean,  et  al.  Since  1990s, 
many  works  with  DNS  have  been  done  in  studying  the  passive  scalars  of  turbulent  flow.  In  recent 
years  some  results  of  DNS  for  passive  scalars  with  relative  high  Schmidt  number  (Sc  >  1)  are  reported 
[5'8l  All  those  DNS  results  of  passive  scalars  are  for  the  incompressible  turbulent  flows,  but  the  scalars 
in  compressible  turbulent  flows  is  more  interested  in  aeronautics  and  astronautics,  for  example,  the 
mixing  of  fuel  and  air  in  supersonic  ramjet  is  a  typical  problem  of  scalars  in  compressible  turbulent 
flow. 

In  this  paper  the  passive  scalars  in  decaying  compressible  turbulence  are  solved  with  direct 
numerical  simulation  by  using  7th  order  accuracy  upwind  difference  scheme  and  8th  order  accuracy 
group  velocity  control  (GVC8)  scheme.  The  start  Reynolds  number  Re=  72,  the  turbulent  Mach 
numbers  Ma=0.2-0.9,  and  the  Schmidt  numbers  of  passive  scalars  Sc=2-10  are  used  in  computation.  In 
order  to  validate  the  computed  results,  the  numerical  experiments  are  made  with  different  simulation 
methods  but  with  same  fluid  parameters.  In  order  to  check  if  the  small  structures  we  are  interested  in 
are  captured,  the  same  problems  are  solved  with  mesh  doubling.  Numerical  experiments  show  that  the 

results  given  in  this  paper  are  reliable.  The  Batchelor  k  1  range  in  scalar  spectrum  is  found  in  our 
simulations.  The  effect  of  compressibility  on  the  flow  structures  is  discussed.  The  results  show  that  the 
high  wave  number  spectrum  decays  faster  with  increasing  the  turbulent  Mach  number.  From  our 
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computed  results  it  is  found  that  in  the  compressible  turbulence  there  is  also  the  extended  self¬ 
similarity  (ESS)  for  the  passive  scalars. 

The  methods  used  in  the  computation  are  presented  in  section  2,  the  validation  of  results  is  given 
in  section  3,  and  the  computed  results  are  analyzed  in  section  4. 


2  Numerical  method 

The  compressible  Navier-Stocks  equations  are  written  in  the  vector  form  as  follows 

0U  8E  8F  8G  8EV  8¥v  8GV 

8t  8x  8y  8z  8x  8y  8z 

and  the  transport  equation  for  the  passive  scalar  expressed  as 

^+d^+d^+d^  =  _^(^%)  +  ^dg^)+d_{judg))  (2) 

dt  dx  3y  dz  Re  Sc  dx  dx  3y  3y  dz  dz 

where  u  =  [p,pu,pv,pw,e]T  ,  E,F,G,EV,FV,GV  are  the  convection  terms  and  the  viscous  terms[9], 
respectively.  The  viscosity  coefficient  jj  is  obtained  from  the  Sutherland  equation,  g  is  the  passive 
function,  and  Sc  is  the  Schmidt  number. 

Periodic  boundary  conditions  and  uniform  meshes  are  used  in  x-,  y-  and  z-  directions.  The  initial 
conditions  of  the  velocity  components,  the  pressure  and  the  temperature  are  given  as  in  Ref.  [3],  and 
the  initial  condition  for  the  passive  scalar  is  given  in  the  same  way  as  that  of  u  component  of  the 
velocity  vector. 

The  viscous  terms  of  the  Navier-Stokes  equations  are  discretized  with  8th  order  center  difference. 
The  flux  vector  splitting  is  adopted  for  the  convection  terms.  The  7th  order  accurate  upwind  difference 
approximation  in  Ref.[3]  and  the  GVC8  in  Ref.[4]  are  used  to  discretized  the  convection  terms. 

Consider  the  approximation  for  the  first  derivative  in  the  convection  terms 

Uj  =  C^/+l/2  —  I  ^ 

For  the  positive  flux  the  following  backward  difference  approximation  is  used 

9 

Fj+ 1/2  =  ^  bkUj+5-k 

k=\ 


|  when 

A u  j~\n 

< 

AWy+l/2  | 

\bk~)  when 

An  j-l/2 

> 

All  /+1/2 

where  A u-  =Uj+i  ~uj-\ 


The  coefficients  b[l)  and  b are  given  in  the  following  table. 


Table  1.  Values  of  b[l)  and  b[2) 


k 

l 

2 

3 

4  5 

6 

7 

8 

9 

h™ 

17 

283 

53 

6269  4429 

10531 

4253 
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3 

7000  ’ 

21000  ’ 

21000  ’ 

21000  ’  4200  ’ 

21000’ 

21000  ’ 

7000  ’ 

500 

b[2) 

4 

893 

4063 

14501  2371 

2299 

137 

31 

1 

875  ’ 

21000  ’ 

21000’  21000  ’  4200  ’ 

21000  ’ 

21000  ' 

7000  ’ 

1000 

For  the  negative  flux  the  following  forward  difference  approximation  is  used: 


9 

Fj+ 1/2  akUj+A-k 

k= 1 


where  ak 


b\o-k  when 

AMy-i/2 

< 

AMy+l/2 

bio-k  when 

Awy-_i/2 

> 

AWy+1/2 

(A:  =  1,2-9) 
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The  difference  approximation  for  the  passive  scalar  is  the  same  as  that  for  the  Navier-Stokes 
equations.  To  save  computing  resource  the  Reynolds  number  in  computation  is  not  very  high,  and  the 
grid  number  for  the  passive  scalar  is  doubled  in  each  direction  comparing  with  the  grid  number  for 
other  flow  parameters.  8th  order  Langrage  interpolation  is  used  to  transport  the  value  of  coarse  mesh  to 
the  fine  mesh. 

The  parameters  used  in  the  DNS  cases  are  shown  in  the  following  table 


Table  2.  Computing  parameters  of  each  DNS  cases 


CASE 

Re  a 

Mt 

Sc 

Numerical 

Method 

Grid  number  for 
fluid 

Grid  number  for 
passive  scalar 

D1 

72 

0.2 

5 

UD7 

256X256X256 

512X512X512 

D2 

72 

0.5 

5 

UD7 

256X256X256 

512X512X512 

D3 

72 

0.7 

5 

UD7 

256X256X256 

512X512X512 

D4 

72 

0.9 

5 

GVC8 

256X256X256 

512X512X512 

El 

72 

0.5 

2 

UD7 

256X256X256 

512X512X512 

E2 

72 

0.5 

10 

UD7 

256X256X256 

512X512X512 

E1T 

72 

0.5 

2 

UD7 

128X128X128 

256X512X512 

D2Ta 

72 

0.5 

5 

UD7 

128X128X128 

256X256X256 

D2Tb 

72 

0.5 

5 

GVC8 

256X256X256 

512X512X512 

In  Table  2  ReA  is  the  initial  Reynolds  number,  Mt  is  the  initial  turbulent  Mach  number,  and  Sc  is  the 

Schmidt  number.  UD7  in  Table  2  denotes  7th  order  upwind  difference  method  [3]  and  GVC8  denotes  8th 
order  group  velocity  control  scheme. 

The  initial  Reynolds  number  equal  to  72  for  all  cases  in  computation,  the  Schmidt  number  is  in 
the  range  2-10,  and  the  initial  turbulent  Mach  number  is  in  the  range  0.2-0. 9.  The  flow  parameters  for 
the  case  E1T  in  Table  2  are  the  same  as  for  the  case  El,  and  the  parameters  for  D2Ta,  D2Tb  are  the 
same  as  for  D2.  These  three  cases  E1T,  D2Ta  and  D2Tb  are  used  to  validate  numerical  results. 

Programs  are  coded  by  using  MPI  Fortran  77,  and  the  DNS  results  are  computed  with  the  LSSC2 
computer  of  LSEC  (State  Key  Laboratory  of  Scientific  and  Engineering  Computing).  Averaged 
performance  are:  11  seconds/time-step  for  case  E1T  and  D2Ta  (8  CPUs  in  8  nodes),  18.6 
seconds/time-step  for  cases  D1-D3  and  E1-E2  (64  CPUs  in  32  nodes),  20.5  seconds/time-step  for  case 
D4  and  D2Tb  (64  CPUs  in  32  nodes). 

3.  Validation  of  computed  results 
3.1  The  flow  parameters 

All  the  flow  parameters  obtained  from  numerical  simulations  have  been  validated  by  the  authors  in 
papers  [3-4], 

In  Fig.l  is  shown  the  time  history  of  normalized  average  kinetic  energy  K(t)/  K0  and  velocity 
derivative  skewness  for  Case  D2.  Horizontal  ordinate  in  Fig.l  denotes  the  normalized  time  tlr . 

Where  r  =  J—  (2n)x/4  k^7'2  is  the  large-eddy-tumover  time  at  t=0  [2l  The  symbols  in  this  figure  are 

V  A 

the  results  given  in  Ref.  [2]  by  using  10th  order  accurate  Pade  schemes  with  the  same  computing 
parameters,  including  the  initial  Reynolds  number  and  initial  turbulent  Mach  number,  as  these  in  Case 
D2.  Fig.l  shows  that  the  results  of  this  paper  agree  well  with  the  results  of  Ref.  [2],  and  this  validated 
our  results  of  flow  fields. 

The  numerical  experiments  for  the  case  with  Reynolds  number  Re/l=72  show  that  the  grid 
number  128X128X128  is  enough  for  obtaining  the  accurate  results.  In  this  paper,  the  gird  systems  128 
X128X128  and  256X256X256  are  used  for  simulating  the  flow  field  parameters. 
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Figure  1.  Time  history  of  normalized  kinetic  energy  and  Figure  2.  Time  history  of  RMS  passive  scalars  for  case  El 
velocity  derivative  skewness  for  case  D2  and  the  results  of  (512*512*512)  and  E1T  (256*256*256) 

Samtaney  et  al. 

3.2  The  passive  scalars 

We  will  compare  the  statistics  of  passive  scalars  in  each  DNS  cases  by  using  different  numerical 
methods  and  with  different  grid  sizes.  Numerical  experiments  show  that  the  statistics  is  not  sensitive 
with  different  numerical  methods  or  grid  size  if  the  grid  number  is  large  enough  (more  than  or  equal  to 
128X128  X128).  This  means  that  the  small  structures  we  are  interested  in  are  captured. 

Figure  2  shows  the  time  history  of  RMS  passive  scalars  for  Case  El  and  E1T.  The  computing 
parameters  and  numerical  method  are  the  same  for  these  two  cases,  but  the  grid  number  are  512X512 
X512  and  256X256X256  for  these  two  case  respectively.  The  RMS  of  passive  scalar  is  defined 

as:  Rg(t)  =  grms(t) / grms(0) ,  grms  =^<  g2  >  ■  From  fig.2  we  can  see  that  the  two  results  agree 

well  with  each  other.  In  Fig.  3  is  shown  the  time  history  of  RMS  passive  scalar  for  case  D2  and  D2Ta, 
from  which  we  also  can  see  that  grid  size  is  small  enough  for  simulating  physical  phenomena  we  are 
interested  in.  This  experiment  tells  us  that  the  error  due  to  grid  size  is  negligible. 

Fig.4  shows  the  time  history  of  RMS  passive  scalar  for  case  D2  and  D2Tb.  The  flow  parameters 
and  the  grid  numbers  in  computation  for  these  two  cases  are  the  same,  but  the  numerical  methods  are 
different.  The  7th  order  upwind  difference  scheme  is  used  for  the  case  D2,  and  GVC8  scheme  is  used 
for  the  case  D2Tb.  From  this  figure  we  can  see  that  the  results  of  these  two  cases  agree  well  with  each 
other.  This  numerical  experiment  also  validate  our  numerical  results. 

As  it  is  known  that  the  grid  number  in  each  direction  needs  to  be  the  magnitude  of  Re34  where 
Pe-  Rex  Sc  is  the  Pe  number.  We  can  infer  that  the  grid  number  for  the  case  D4  is  enough  because 
the  Pe  number  for  case  D4  is  the  same  as  for  the  case  D1-D3.  We  can  also  infer  that  the  grid  number 
for  E2  is  sufficient  because  the  grid  number  for  E2  is  doubled  comparing  with  the  case  D2Ta  in  each 
direction.  This  grid  number  is  more  than  the  needed. 


Figure  3.  Time  history  of  RMS  passive  scalar  for  case  D2  r  ^ 

J  r  Figure  4.  Time  history  ot  RMS  passive  scalar  tor  case  D2 

(512*512*512)  and  D2Ta  (256*256*256)  (GVC8  scheme)  and  D2Tb  (7th  order  upwind  difference 

scheme) 
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4.  Analysis  of  Numerical  Results 

We  will  study  the  flow  field  and  passive  scalar  field  for  each  DNS  case  at  t/r  =  1.  The 
parameters  of  these  fields  are  showed  in  the  following  table.  Turbulence  after  tlr-  0.8  may  be 
considered  realistic  [2'3].  The  results  obtained  at  t  /  v  =  1  are  analyzed. 


Table  3.  The  parameters  of  the  flows  and  passive  scalar  fields 


Field 

Number 

DNS 

CASE 

tlr 

Rez 

Mt 

Sc 

FD1 

D1 

1.0 

39.4 

0.17 

5 

-0.44 

-0.63 

FD2 

D2 

1.0 

38.5 

0.42 

5 

-0.40 

-0.63 

FD3 

D3 

1.0 

36.3 

0.42 

5 

-0.42 

-0.51 

FD4 

D4 

1.0 

33.4 

0.72 

5 

-0.69 

-0.57 

FE1 

El 

1.0 

38.5 

0.42 

2 

-0.40 

-0.64 

FE2 

E2 

1.0 

38.5 

0.42 

10 

-0.40 

-0.62 

where  Su 


<u  > 

2  3  n 

<  u  > 


ug 


<u2g> 


.  2  .  .  2 

<  U  ><  g 


> 


1/2 


O 


4.1  Spectral  analysis  of  passive  scalar 

Batchelor  (1959)  found  that  there  is  AT1  range  in  the  energy  spectrum  of  passive  scalar  when  the 

Schmidt  number  is  much  large  than  1.  In  later  researches  it  was  found  that  “AT1  range”  in  the 
spectrum  exists  only  when  the  Schmidt  number  large  than  1  (need  not  much  large  than  1).  But  in  some 
experiments  “  k  1  range”  for  the  passive  scalar  spectrum  has  not  been  found  yet  |fi|.  In  all  these 
research  works  the  passive  scalars  are  discussed  only  for  the  incompressible  turbulence,  and  the 
effects  of  compressibility  on  passive  scalars  are  still  unknown. 

The  energy  spectrum  of  passive  scalar  is  defined  as: 

Eg ( k )  =  4 7k 2 (g(kx ,k2,k3)g* (k{  ,k2,k2 


/V  .  /V  ^  .  .  /V 

where  g  is  the  spectrum  of  g  ,  g  is  the  conjugate  of  g  .  The  shell  average  of  spectrum  at  the  range 
\k  - 1  /  2,  k  + 1  /  2)  is  expressed  as 

(f(kuk2,k3))k=^~  £/(*!,  *2,  *3 ) 

^  k  k-U2<k<k+V2 

where  k  =  yjk2  +  k2  +  k 2  . 


Figure  5.  Energy  spectrum  of  passive  scalar  for  FD1-FD4 


Figure  6.  Energy  spectrum  of  passive  scalar  for  FE1,  FD2 
and  FE2 


Fig.5  shows  the  energy  spectrum  of  passive  scalar  E  for  FD1-FD4  (scalar  fields  at  tlr  =  1  of 

cases  D1-D4).  We  can  find  that  the  spectrum  decays  faster  in  the  range  of  high  wave  number  with 
turbulent  Mach  number  increasing.  This  is  because  that  more  kinetic  energy  will  be  dissipated  by 
shocklets  with  turbulent  Mach  number  increasing. 

Fig.6  shows  the  energy  spectrums  of  passive  scalars  for  the  cases  FE1,  FD2  and  FE2.  The 
Schmidt  number  for  these  three  cases  is  2,  5  and  10  respectively.  This  figure  shows  that  the  energy 
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spectrum  of  passive  scalar  decays  slower  with  the  increasing  Schmidt  number.  The  amplitude  of  the 

—8 

highest  wave  number  of  the  spectrum  is  about  10  which  shows  the  sufficient  accuracy  of  our 
simulations. 

Fig.7  shows  the  normalized  energy  spectrum  of  passive  scalars  for  cases  FE1  FD2  and  FE2.  The 
horizontal  ordinate  is  kriB  and  the  vertical  ordinate  is  (kr/B)E g /(%(u / s)1" rfB)  where  r/B  is 
Batcher  scale,  £  and  x  's  the  dissipate  rate  of  kinetic  and  passive  scalar  respectively.  The  normalized 
energy  spectrum  (krjB)E g  /(x(v / £)l/2  rjB)  remains  unchanging  in  a  relative  long  interval.  This 

result  means  that  exists  relation  E  ~  k  1  in  this  interval,  and  this  Batchelor  k  1  spectral  is  universal 
for  all  these  cases  considered. 


Figure  7.  Normalized  energy  spectrum  of  passive  scalar  of 
FE1,  FD2  and  FE2. 


4.2  Coherent  structure  of  passive  scalar  in  compressible  isotropic  turbulence 

The  sheet-like  structures  in  passive  scalar  in  incompressible  turbulence  are  reports  in  some  recent 
research  works,  but  there  still  no  reports  about  coherent  structures  of  passive  scalar  in  compressible 
turbulence. 
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where  0rms  is  the  root-mean-square  of  the  scalar 


Figure  8  and  Figure  9  shows  the  isosurface  of  9  —  2  for  FD2  and  FD4,  respectively.  The 
visualizing  domain  is  the  1/64  of  the  full  computing  domain.  The  two  figures  show  the  sheet  structure 
of  passive  scalar  in  compressible  turbulence. 


Figure  8.  Isosurface  of  normalized  scalar  gradient  of  Figure  9.  Isosurface  of  normalized  scalar  gradient  of 
FD2  FD4 
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4.3  Scaling  analysis  of  passive  scalar  in  compressible  isotropic  turbulence 

Scaling  law  is  a  hot  topic  in  turbulence  physics  in  recent  years.  Extended  self-similarity  (ESS)  is 
one  of  the  important  developments  in  measurement  of  the  scaling  exponents.  She  et  al  developed  S-L 
model  of  scaling  law,  which  has  been  tested  by  several  experiments  and  numerical  dates.  The  scaling 
law  of  the  passive  scalar  in  turbulent  flow  has  been  reported  in  recent  years.  Zhou  et  al|s|  studied  the 
scaling  law  of  passive  scalar  for  the  incompressible  isotropic  turbulence.  But  there  is  still  no  report  for 
scaling  law  of  passive  scalar  for  the  compressible  turbulence. 

Fig.  10  are  given  the  plots  of  p-th  order  passive  scalar  structure  functions  as  function  of  3rd  order 

passive  scalar  structure  function  for  the  case  FD2.  The  horizontal  ordinate  is  log  <|  dgt  |3> ,  and  the 

vertical  ordinate  is  log  <|  5gt  |p>  where  dgl  =  g(x  +  /)  -  g(x)  ,  and  (■)  denotes  the  average  over 

the  all  field.  From  the  plot  we  can  see  clearly  the  linear  variation,  which  means  that  the  power  laws 
between  p-th  order  passive  scalar  structure  function  and  3rd  order  passive  scalar  structure  function. 
This  means  ESS  (Extended  Self-Similarity)  holds  in  the  compressible  passive  scalar  turbulence. 

Fig.l  1  and  Fig.  12  are  shown  the  same  plots  for  the  passive  scalar  field  for  the  case  FD4  and  FE2, 
respectively.  From  these  two  figures  we  can  obtain  the  same  conclusion  as  that  in  figure  10.  The  rates 
of  slope  in  the  plots  in  fig.  10-12  are  the  relative  scalar  exponent  C,  ,  and  the  least  square  approximate 

method  is  used  to  calculate  these  slope  rates. 


log  <i<%,r> 


Figure  10.  Plots  of  p-th  order  passive  structure  functions  as 
functions  of  3rd  passive  structure  in  FD2 


log  <|  Sg,\p> 


as  functions  of  3rd  passive  structure  in  FD4 


Fig.  13  shows  the  relative  scaling  exponents  for  the  cases  FD1-FD4.  The  dash  line  denotes  p/3 
law.  The  exponents  lower  than  p/3  denote  the  strong  intermittence  of  the  passive  scalar  field.  This 
figure  indicates  that  the  compressibility  has  little  effect  on  passive  scalar’s  scaling  exponents,  for 
example,  the  largest  difference  of  <^6  is  not  more  than  4%  and  the  largest  difference  of  (T  is  not  more 
than  5%  among  the  four  cases. 

Fig.  14  shows  the  relative  scaling  exponents  for  FE1,  FD1  and  FE2  with  Schmidt  numbers  2,  5 
and  10  respectively.  This  figure  shows  scaling  exponent  become  smaller  with  increasing  Schmidt 
number.  It  means  that  the  passive  scalar  become  more  intermittent  with  the  increasing  Schmidt 
number. 


Figure  12.  Plots  of  p-th  order  passive  structure  functions  as  functions  of  3rd  passive  structure  in  FE2 
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Figure  14.  Relative  scaling  exponents  of  passive  scalar  for 
FE1,  FD2  and  FE2 

5.  Conclusion 

1.  There  is  clear  Batchelor  k  l  range  in  the  energy  spectrum  of  passive  scalar  in  the  compressible 
turbulent  flow,  and  the  compressibility  effect  leads  to  faster  decay  of  energy  spectrum  of  passive 
scalar  in  high  wave  number  ranges. 

2.  The  sheet-like  structures  of  the  passive  scalar  gradient  are  founded  in  compressible  isotropic 
turbulence. 

3.  The  ESS  holds  in  the  passive  scalar  of  compressible  turbulence,  and  the  compressibility  has 
little  effects  on  the  relative  scalar  exponents  of  passive  scalar. 

4.  The  exponents  of  passive  scalar  became  more  intermittent  with  Schmidt  number  increasing. 
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